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Abstract 



We present the results of extensive molecular dynamics computer simulations 

in which the high frequency dynamics of silica, v > 0.5 THz, is investigated 

in the viscous liquid state as well as in the glass state. We characterize the 

. ■ properties of high frequency sound modes by analyzing Ji{q,v) and Jt(g, z/), 

C^ . the longitudinal and transverse current correlation function, respectively. For 

^ I wave-vectors q > 0.4 A^^ the spectra are sitting on top of a flat background 

i-rt I which is due to multiphonon excitations. In the acoustic frequency band, 

^ ■ i.e. for 1/ < 20 THz, the intensity of Ji{q,v) and Jt{q,i') in the liquid and 

the glass approximately proportional to temperature, in agreement with the 

harmonic approximation. In contrast to this, strong deviations from a linear 

scaling are found for ly > 20 THz. The dynamic structure factor S{q, v) 

^ \ exhibits for q > 0.23 A^^ a boson peak which is located nearly independent 

Ti^lj- ■ of q around 1.7 THz. We show that the low frequency part of the boson 

^T ■ peak is mainly due to the elastic scattering of transverse acoustic modes with 

frequencies around 1 THz. The strength of this scattering depends on q and 

0^ i is largest around q = 1.7 A~^, the location of the first sharp diffraction peak 

- ' in the static structure factor. By studying S{q, u) for different system sizes we 

C^ ' show that strong finite size effects are present in the low frequency part of the 

boson peak in that for small systems part of its intensity is missing. We discuss 

i-A . the consequences of these finite size effects for the structural relaxation. 
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I. INTRODUCTION 

The investigation of the vibrational dynamics of supercooled liquids and glasses is a 
challenging task since these systems do not have the property of translational invariance as 
it is the case for crystals. Of special interest is the region of intermediate wave-vectors at 
which collective excitations, i.e. longitudinal and transverse sound waves, begin to experience 
strongly the structural disorder. Molecular dynamics (MD) simulations are well suited to 
study vibrational features at these intermediate wave-vectors, say with magnitude q > 
0.1 A~^, corresponding to frequencies u in the THz band. This paper is concerned with the 
simulation of the vibrational dynamics of amorphous silica, which is the prototype of a so- 
called strong glassformer |jl[. Its structure exhibits a medium-range order in that it forms 
a disordered network of Si04-tetrahedra leading to a peak in the static structure factor 
around q = 1.6 A~^ 0. In recent years the high frequency dynamics of silica has been the 
subject of an intense debate because its Raman and neutron scattering spectra ||^,|] show 
a so-called boson peak around 1 THz, which is also present in many other glassformers 
but normally gives a less intense contribution to the spectra than in silica. This feature 
appears also in experiments as an excess over the Debye density of states or, equivalently, 
over Debye's T^-law in the specific heat around the temperature T = 10 K [§,H. 

Recently it has been shown by means of an inelastic X-ray scattering experiment by 
Benassi et al. that in sihca propagating longitudinal sound modes persist up to 0.35 A~^, 
which corresponds to frequencies well above the location of the boson peak. Therefore Be- 
nassi et al. argued that the boson peak has its origin in these propagating sound modes. In 
contrast to this suggestion, Vacher et al. [pHlOU found evidence from the available spectro- 
scopic data that the boson peak is due to the strong scattering of acoustic modes by the 
disorder thus regarding the boson peak as a manifestation of these strongly scattered modes. 

Simple models have been suggested to explain the origin of the boson peak. From the soft 
potential model [lTIl , p!2| the idea is put forward that anharmonic localized vibrations coexist 



with propagating high frequency sound modes in the frequency range around the location 
of the boson peak. In the case of silica these anharmonic soft modes have been related to 
coupled Si04-tetrahedra librations 0]. Wischnewski et al. [|T3| have analyzed their neutron 
scattering data of silica within the soft potential model, and have concluded that the sound 
waves are indeed scattered from such local vibrational modes below 1 THz, whereas above 
this frequency static Rayleigh scattering from the atomic disorder takes place. Schirmacher 
et al. 0] have studied a system of coupled harmonic oscillators with a random distribution 
of force constants. In this model they have found an excess over the Debye behavior in the 
density of states which they have interpreted as an analogon to the boson peak feature in 
real structural glasses. In agreement with this model, Sokolov []T3| proposed that the boson 
peak is related to the strong scattering of acoustic like vibrations by fluctuations of the 
elastic constants. 

A feature which shares many properties with the boson peak is also found within the 
mode-coupling theory (MCT) of the glass transition: In the ideal glass state where all 
particles are trapped in the cages formed by their neighbors, the spectrum of the density- 
density correlation function is a superposition of harmonic oscillator spectra which is due 
to the variety of cages in which the particles are trapped |]16|. It is remarkable in this 



context that the whole light scattering spectra of glycerol, including the boson peak, have 



been successfully described within a schematic MCT model JT^, and very recently Gotze 
and Mayr have shown that deep inside the glass state, i.e. at temperatures well below the 
MCT temperature T^, the theory predicts dynamical features which are very reminiscent to 
the boson peak 



In the last three years MD simulations tried to give insight into the vibrational dynamics 



of silica [pJH 2^ , and other network forming glasses like ZnCl2 [^ • Most of these investiga- 
tions have analyzed the dynamics within the harmonic approximation, i.e. by determining 
the eigenvalues and eigenvectors from the diagonalized dynamical matrix. Although the full 
information of the vibrational part of the dynamics is given by the eigenmodes (of course 
only within the harmonic approximation) the origin of the boson peak remains a puzzle. 
One reason for this is the smallness of the system sizes (20-40 A) which have been used 
in the forementioned studies, which has the effect that parts of the boson peak are missing 
(see below). A second reason is the difficulty of analyzing the boson peak feature in terms 
of eigenmodes since, as we will discuss in detail below, in the case of silica the couphng of 
modes with different q is essential for this feature. 

In order to avoid these problems we use in the present work a large system size and 
calculate the exact current and density correlation functions in order to investigate their 
dependence on wave-vector q and frequency v. Therefore, we are able to study the tem- 
perature dependence of the high frequency dynamics of silica in the liquid state as well as 
in the glass state. Moreover, we are able to give insight into the relationship between the 
vibrational dynamics and structural relaxation in the silica melt. The rest of the paper is 
organized as follows: In the next section we give an overview of the main computational 
details. In Sec. |T1| we discuss the vibrational dynamics of our silica model by means of the 
current and density correlation functions. In Sec. |I^ we summarize and discuss the results. 

II. MODEL AND DETAILS OF THE SIMULATIONS 

The model potential we use to describe the interactions between the ions in silica is the 
one proposed by van Beest, Kramer, and van Santen (BKS) [^ which has the following 
functional form: 

0(r) = ^^ + A,/5exp(-5,^r)-^ «,/?e[Si,0]. (1) 

Here r is the distance between an ion of type a and an ion of type (3. The values of the 
parameters Aq,^, i?^^ and Cap can be found in the original publication. The Coulombic part 
of the potential was evaluated by means of Ewald sums for which further details can be found 
elsewhere [0. In recent simulations [^-31| it has been shown that the BKS potential (|l]) 



reproduces many static and dynamic properties of real silica very well and thus it can be 
considered as a reliable model for this material. 

We have simulated a system with 8016 ions. The size of the simulation box was fixed to 
48.365 A corresponding to a density of 2.37 g/cm^. Thus, the smallest wave-vector of our 
simulation has magnitude g = 0.13 A~^. In order to study finite size effects we have done 
also simulations for smaller systems, and the details of these simulations are given below. In 
the following we will investigate the fully equilibrated liquid state at T = 6100 K, 3760 K, 
and 2750 K and the glass state at T = 1670 K, 1050 K, and 300 K. In the liquid state we 



equilibrated the system first at each temperature in the NVT ensemble at each temperature, 
and after that we started microcanonical simulations by means of the velocity form of the 
Verlet algorithm. During the equilibration the temperature was kept constant by using a 
stochastic collision algorithm. The time step we used was 1.6 fs, and in order to improve 
the statistics we simulated at each temperature two independent runs. At T = 2750 K the 
length of the equilibration runs was 13 million time steps followed by the microcanonical 
production runs over 12 million time steps, which corresponds to a real time of 20 ns. 
During the two production runs we have stored on a linear time scale 30 configurations each 
which have subsequently been used as the starting configuration of a new simulation for 
investigating the high frequency dynamics. We mention that the pressure at T = 2750 K is 
around 0.9 GPa. Further details on the simulation of the liquid state can be found elsewhere 
27| . The starting-point for producing the glass state were two equilibrated configurations 



at T = 2900 K at which the equilibration time was 4 million time steps (6.5 ns real time). By 
coupling the system to an external heat bath the temperature was then decreased linearly 
in time within one million time steps to K. This corresponds to a cooling rate of about 
1.8 ■ 10^^ K/s. During the cooling procedure we stored configurations at the temperatures 
mentioned above which we used as starting configurations in order to anneal the system for 
5 ■ 10^ time steps at constant temperature. Afterwards we propagated the system over 5 ■ 10^ 
time steps in the microcanonical ensemble and stored configurations every 10^ time steps. 
Thus at the end we had at each of the three temperatures in the glass state 22 starting 
configurations for the investigation of the high frequency dynamics. The pressure for our 
glass structures is 0.52 GPa at T = 300 K, 0.69 GPa at T = 1050 K, and 0.8 GPa at 
T = 1670 K. 

In this paper we are mainly interested in frequency dependent correlation functions. 
Therefore time Fourier transformations have to be calculated which we have done by means 
of the Wiener-Khinchin theorem. It says that the Fourier transformation of a correlation 
function C(t) = {x(t)x{0)) {x{t): density, longitudinal current, transverse current) is given 
by the power spectrum 2'(z/) = \a{i')\ where a(z/) denotes the Fourier transform of the 
time series x{t). The time series were transformed via fast Fourier transformation whereby 



we applied a Welch window function |^. Usually we have calculated the time series for 
the density and the currents over 8192 time steps (13.4 ps real time) by using the afore- 
mentioned starting configurations. This results in a frequency resolution of about 0.1 THz. 
The reliability of the Fourier transformation was tested by calculating also time series over 
16384 time steps and in these test cases we have found indeed identical spectra, at least for 
u > 0.3 THz. 



III. RESULTS 

A. Current correlations 

In this section we analyze the vibrational features of our silica model by means of the 
longitudinal and transverse current correlation function Ji^q,^) and Jt^q,!^), respectively, 
which depend on the magnitude of the wave-vector q and the frequency u. These are 
defined as |Q 



1 f°° 
Jaiq,J^) = j^J c?t exp(i27rz/t) (j<^(q,t) ■j„(-q,0)) (2) 

where the longitudinal part {a = I) and the transverse part {a = t) of the total current 



K^^t) =J2rk{t) exp {iq-Vk (t) ) (3) 

k 



are given by 



.. ,s qq-j(q,^) /.^ 

3m, t) = — , (4) 

. . ,^ V +\ qqj(q^^) /.^ 

jt q, ^ = J q, t) ^ — . (5) 

Fig. shows Ji{q, u) and Jf (g, u) for different values of q up to 1.0 A~^ at the temperature 
T = 2750 K. We have plotted only the functions for the oxygen-oxygen correlations because 
the silicon-silicon and the silicon-oxygen correlations exhibit qualitatively the same behav- 
ior, which is reasonable for such small wave-vectors. Note that even at the relatively high 
temperature T = 2750 K, our Si02 model is quite viscous having a viscosity of about 380 P, 
and moreover, that this temperature is well below the critical temperature of mode-coupling 



theory, which is at 3330 K [£7|. Thus within the framework of MCT we are indeed probing 
the system deep in the glass regime. In Fig. 0a we show J/(g, z/) and Jj(g, u) in the frequency 
range between 0.4 and 1.6 THz for the four lowest q values of our simulation, g = 0.13 A~^, 
0.18 A~^, 0.23 A~^, and 0.26 A~^. At g = 0.13 A~^ we recognize that there are two peaks, 
corresponding to the longitudinal and the transverse part of the current, which are well sep- 
arated from each other. For increasing wave-vectors these peaks move to higher frequencies 
whereby their width becomes so large that they overlap more and more with each other. In 
the following we call the excitations corresponding to these peaks high frequency longitudi- 
nal acoustic (LA) modes and high frequency transverse acoustic (TA) modes, respectively. 
From the figure we see that the TA excitations give the most important contribution to the 
current spectra in that their amplitude is about a factor 6-8 higher than that of the LA 
excitations. In the wave-vector range in which the LA and TA modes hybridize one would 
expect that plane waves are no longer eigenmodes, and in the simulation of Taraskin and 
Elliott it has indeed been shown explicitly that a longitudinal or transverse plane wave with 
a q value around 0.2 A~^ decays into a final state which can be characterized as a super- 
position of plane waves with different wave-vectors and polarizations, but with the same 
frequency ppj] . 

At g = 1.0 A^^ (Fig. |l|b) the current correlation functions are qualitatively different from 
those discussed so far at lower g: In the transverse part one observes a plateau between 3 
and 11 THz, and in the longitudinal part the LA peak around 16 THz seems to be sitting 
on top of a flat background. In order to describe in more detail the change in the shape of 
the spectra that occurs at intermediate values of g, we show in Fig. ^ Ji{q, v) and Jt(g, u) 
for the 0-0, Si-0, and Si-Si correlations for g up to 1.7 A"^. At g = 0.47 A~^ we observe 
in J/(g, v) for the 0-0 correlations, apart from the LA peak around 6 THz, a peak around 
z/ = 26 THz corresponding to an optical excitation. Moreover, the intensity of the whole 
spectrum seems to be enhanced in that the LA and the optical peak sit on top of a flat 



background. If q is increased to 1.4 A~^ the LA peak moves to larger frequencies whereby a 
shoulder around v = 2 THz gets more and more pronounced. Also at g = 1.7 A~^ there is a 
LA peak but now its position has moved back to i/ = 17 THz, i.e. a smaller frequency than 
at g = 1.4 A~^. In the case of the Si-0 correlations (Fig. ^) the essential difference to the 
0-0 correlations is the negative amplitude of the LA peak for g > 1.4 A~^ which indicates 
an antiphase motion of the silicon and oxygen atoms. The curves for the Si-Si correlations 
(Fig. He) show essentially only one difference compared to those for the 0-0 correlations in 
that the optical band has a higher weight in the spectrum than the LA excitations. This is 
due to the fact that the silicon atoms are bonded stronger in the tetrahedral network than 
the oxygen atoms, and thus on small length scales more localized motions have a higher 
weight in the case of the silicon atoms which corresponds to frequencies in the optical band. 
Also in the transverse case for the 0-0 correlations (Fig. Qd) the whole spectrum sits 
on top of a flat background. The intensity of the TA peak around 3 THz decreases with 
increasing q whereas there is an increase in the intensity around 9 THz. As a result a 
broad fiat band is obtained for z/ < 17 THz. In contrast to the 0-0 correlations, Jt{q-,v) 
for the Si-0 correlations (Fig. |^e) shows a strong overall decrease of the intensity if q is 
increased from 1.0 A~^ to 1.7 A^^. This can be easily understood because at g = 1.7 A~^ 
the current correlation functions measure to a great extent the single particle motion, and 
therefore the relative motion of the silicon and oxygen atoms gives only a small contribution 
to the spectra. The most remarkable feature in Jj(g, z/) for the Si-Si correlations is again 
that, compared to the 0-0 correlations, the optical excitation around 20 THz has a larger 
amplitude than those of the acoustic band for g > 1.0 A~^. The essential result which is 
shown in the Fig. |^ is that for intermediate values of g the whole spectrum is placed on 
top of a fiat background. A similar feature has also been found by Mazzacurati et al. 



in a Lennard- Jones system. These authors have identified the fiat background directly 
in the spectra and in the participation ratio which measures the number of particles that 
contribute to the eigenmodes at a certain frequency. At the low frequency edge of the 
density of states the participation ratio has values expected for localized modes. Such a 
behavior of the participation ratio has also been found in the case of silica ||3^ . Mazzacurati 
et al. have explained this behavior by showing that the eigenvectors for low frequencies can 
be represented by a few long-wavelength standing waves plus a random contribution where 
the random contribution is seen in the spectrum as the fiat background. In a phonon picture 
one can interpret the fiat background as the contribution of multiphonon excitations. 

By reading off the peak maxima [^ in Ji{q,i>) and Jt{q,v) corresponding to the lon- 



gitudinal and transverse acoustic modes one gets dispersion like branches z/;(g) and t't(g) 
which are shown in Fig. |^ for T = 2750 K and in Fig. ^d for T = 300 K. It is remark- 
able that vi{q) and vt{q) exhibit essentially the same behavior in the viscous liquid state 
(T = 2750 K) and the glass state (T = 300 K). This shows that in this (high) frequency 
window there is no relevant difference between a viscous liquid and a glass which gives sup- 



port to the idea of Ref. |0 that in this frequency range the viscous liquid can be treated 

like a glass. Furthermore, we note that both functions look very similar as in simple liquids 

5B[ : The longitudinal branch //^(g) has a periodic structure with a minimum located around 



gm = 2.8 A~^, which is the location of the second sharp diffraction peak in the static struc- 
ture factor and which corresponds to length scales of intratetrahedral distances [^ . Thus, 
gni/2 can be interpreted, in analogy to crystals, as a quasi Brillouin zone. The minimum in 



z/;(g) at gm can be easily understood since the particles tend to favor relative separations of 
2vr/gin, and therefore, at these wavelengths it costs a relatively small amount of energy to 
excite a collective mode corresponding to a relatively small frequency. That the minimum 
in i'i{q) is not observed at g = 1.7 A~^, the location of the first sharp diffraction peak in the 



static structure factor ^^ , is due to the fact that this q value corresponds to length scales of 
connected Si04-tetrahedra, a structural unit which is less stiff than one tetrahedron itself. 
The behavior of i^iiq) is in agreement with the findings in a neutron scattering experiment 
by Aral et al. ||3^, and was also found in the computer simulations of Taraskin and Elliott 
| |35| |. The transverse branch z/((g) becomes rather flat for q > 0.9 A~^ which is an indication 
of the overdamped character of the TA excitations at these wave-vectors. 

Also included in Figs. |a and |]b are fits of the form //^(g) = Caq/i^n), where q and q 
denote the longitudinal and the transverse high frequency sound velocity, respectively. The 
values for c^ obtained from these fits are given in the figures. We recognize that for q up 
to around 0.4 A~^ this linear dispersion law holds, which is expected for propagating sound 
waves at sufficiently small q. We have determined the longitudinal and transverse sound 
velocity for all temperatures considered by calculating Ca = 27rz/a/g for the two lowest q 
values of our simulation g = 0.13 A~^ and g = 0.18 A~^. The sound velocities obtained in 
this way are shown in Fig. ^ as a function of temperature. Note that c^ as determined from 
g = 0.13 A~^ and from g = 0.18 A~^ differ by less than 7 % from each other, which shows 
that these wave-vectors are small enough to determine c^ reliably. From 3760 K to 6100 K 
the longitudinal sound velocity decreases by about 50 %. No data is shown for q at 6100 K in 
the figure because at this temperature only a peak at z/ = is observed. This behavior is in 
agreement with hydrodynamics which predicts that transverse fiuctuations are transported 
diffusively and therefore contribute to the spectrum only with a peak at z/ = 0. We have 
found, however, that even at T = 6100K the restoring forces between the particles are large 
enough to allow the propagation of TA modes for g > 0.35A~^, which can be inferred by the 
observation of a crossover from a peak around z/ = to a peak at finite frequencies in this 
region of g. Also included in the figure are the experimental sound velocities measured by Vo- 



Tanh et al. [Q which are multiplied with the factor J2.2/2.37. This factor takes into account 
that the density of our simulation, 2.37 g/cm^, is slightly different from the experimental 
one, which is 2.2 g/cm'^. With this "correction" the simulation reproduces the experimental 
data very well, both for the longitudinal and the transverse sound velocities. Note that 
the experimental data have been obtained by Vo-Tanh et al. by means of light scattering 
experiments for values of g of the order 10~^ A~^, i.e. about two orders of magnitude below 
the g values of our study. Since, however, it has been shown by Benassi et al. that at 
least the longitudinal sound velocities do not change in this g range, i.e. essentially the same 
value for q is measured in Brillouin scattering experiments and in X-ray scattering up to 
g ~ 0.35 A~^, it is reasonable to compare the values of Cq, from Ref. |^ with our data. 



In order to determine, independent from a model, the width of the peaks corresponding to 
the LA and TA excitations we have plotted the ratio Ja(g, z^)/Jo,max(g, ^) versus z/ — z^Q,(g) 
where Ja,max denotes the amplitude of Jq, (g,z^) at the maximum frequency z/Q,(g). The 
resulting curves are shown in Fig. ^ in the longitudinal case for g values up to 2.0 A^^ and 
in Fig. ^3 in the transverse case for g values up to 0.5 A~^ from which we have read off the 
full width at half maximum rQ,(g). We recognize from Fig. ^ that the broadening in the 
curves for the longitudinal part is non-monotonous for g > 1.0 A~^ whereas the curves in 



the transverse case broaden monotonously, as can be seen in Fig. Pd. 

The hnewidths Ta{q) obtained in this way are shown in Fig. |]for the LA and TA modes. 
We recognize from this figure that a quadratic fit describes the longitudinal half width well 
in the q interval 0.18 A~^ < Q' < 0.5 A~^. Such a behavior has also been found in the 
experiment by Benassi et al. [|^ which was done at T = 1050 K. The linear dispersion for 
i'i{q) and the quadratic law Ti{q) oc q^ support the picture that for q < 0.4 A~^ the system 
behaves like an isotropic elastic medium with respect to the propagation of the bare LA 
phonons. For q > 0.6 A~^ ^i{q) becomes rather fiat up to g = 1.1 A~^. Then this function 
decreases significantly and reaches a minimum in the vicinity of the first sharp diffraction 
peak at g = 1.6 A^^. This is probably due to the fact that the strong spatial correlations 
on the length scale of two connected Si04 tetrahedra decreases the damping of the LA 
excitation. It is interesting that in the q range 0.6 A~^ < q < 2.0 A~^ the width Ti(q) is 
significantly smaller than z/;(g) (see Fig. |]) which means that the LA excitations cannot be 
characterized by a loffe-Regel crossover. Note that no data points are available between 1.1 
and 1.4 A~^ because the LA peak overlaps with the optical band in this q range and thus 
they cannot be identified uniquely. 

From Fig. |^ we also see that the transverse peak width Tt[q) can be described by an 
effective power law with exponent 2.5, Tt{q) oc g^'^. Thus, the TA phonons seems to be 
damped stronger than expected from an isotropic elastic medium which would give an expo- 
nent 2. In the q region above 0.5 A~^ the width Tt{q) becomes larger than the corresponding 
location of the maximum of the peak z/i(g). Therefore, we observe a loffe-Regel crossover in 
the transverse case where the TA excitations lose their propagative character and become 
strongly overdamped. Note that this is the q region for which t't(g) becomes more or less 
flat, in contrast to z/;(g). 

If the current correlation functions would behave as expected from the harmonic ap- 
proximation they would simply scale with temperature within the classical treatment of our 
MD simulation. In order to see to what extent a harmonic approximation holds, we have 
plotted in Fig. J;(g,z/)/T and Ji(g,i/)/T at g = 0.26 A"^ for the 0-0, Si-0, and Si-Si 
correlations. We recognize from these figures that the different peaks corresponding to the 
LA and TA excitations fall very well onto one master curve in the whole temperature range 
3760 K > T > 300 K. Even for T = 3760 K only weak deviations from the curve for 
T = 300 K are visible at low frequencies. In contrast to these acoustic excitations the op- 
tical band, i.e. the excitations for i/ > 20 THz, exhibits strong deviations from a harmonic 
behavior. In the following we denote the different excitations in the optical band as LOl 
and L02 in the longitudinal case and as TOl, T02, and T03 in the transverse case. LOl 
and L02 correspond to the frequencies around 25.5 THz and 35.5 THz, respectively, and 
TOl, T02, and T03 correspond to frequencies around 19 THz, 22.3 THz, and 32.2 THz, 
respectively. (Note that we have been able to identify the different optical branches as a 
function of g by observing different peaks at approximately the same frequency for different 
values of g.) It has been shown by computer simulations 0,^] and experiments 0] that 



the L02 and T02 modes correspond to intra-tetrahedral stretching modes, whereas the op- 
tical modes with lower frequencies are mainly due to bending and rocking motions of larger 
structural units. With decreasing temperature the L02 and T02 peaks shift respectively 
from 32 to 36 THz and from 28.5 to 32 THz and their amplitude increases. The mixing 
of the T02 and L02 modes becomes evident in that a shoulder appears at T = 300 K in 



Ji{q,u) around the frequency of the T02 peak in Jt^q,^) and also in Jti^q^v) a feature is 
seen at the frequency of the L02 mode. From the figure we also recognize that the optical 
excitations at lower frequencies than L02 and T02 have a weaker temperature dependence. 
This behavior is reasonable since these modes are, in contrast to L02 and T02, more ex- 
tended since they originate from the inter-tetrahedral motion. Nevertheless, also the LOl 
peak increases in amplitude with decreasing temperature and shifts to higher frequencies. 
The TOl peak evolves into a double peak structure (TOl and T03) in lowering the temper- 
ature. In Fig. P we show the current correlation functions scaled with temperature at the 
higher wave-vector q = 1.7 A~^. The current correlation functions still scale approximately 
linear with temperature for z/ < 20 THz. For these frequencies the strongest deviations from 
such a behavior are found in the peak corresponding to the LA excitations around 18 THz. 
Moreover the mixing of the L02 excitations with the T02 excitations is much stronger at 
q = 1.7 A~^ than at g = 0.26 A~^. Therefore, at q values around 1.7 A~^ it makes no sense 
to distinguish these excitations as transverse and longitudinal ones. 

Figures ^ show the locations of all the peaks in Ji^til) that correspond to the optical 
excitations that can be identified from the current spectra for the Si-Si and the 0-0 cor- 
relations at the temperatures T = 300 K and T = 2750 K, respectively. Whereas it is 
possible to distinguish at T = 300 K the L02 from the T02 peaks for the whole q range 
0.13 A-i < g < 4.5 A-\ this is not possible at T = 2750 K for g > 2.0 A-^ because a 
strong mixing between L02 and T02 occurs for these values of q. The same is the case at 
T = 300 K but at this temperature the L02 and T02 branches can be distinguished from 
each other in that they form a double peak stucture in the longitudinal and the transverse 
part. In Fig. ^a the mixing contributions at T = 300 K are plotted as the dashed and 
solid lines for the Si-Si and 0-0 correlations, respectively: They show the occurrence of a 
longitudinal mode in the transverse spectrum and, vice versa, the occurrence of a transverse 
mode in the longitudinal spectrum. In this context the denotation of the modes as longitu- 
dinal and transverse ones means that at very low values of q they are expected to be purely 
longitudinal and transverse, respectively. 

For T = 300 K the optical branches at lower frequencies, i.e. LOl, TOl, and T03, can 
be clearly identified in the Si-Si correlations and the 0-0 correlations for q < 0.8 A~^. At 
larger values of q LOl and T03 appear as the dominant contribution in the Si-Si correlations, 
whereas for the 0-0 correlations the spectra are dominated by the T03 excitations in Ji 
and the LOl excitations in Jj. At T = 2750 K the spectra are less pronounced, and we 
can identify only the TOl and the LOl branch apart from the high frequency band above 
30 THz. They are visible in the 0-0 and the Si-Si correlations for q < 0.8 A^^. At higher 
q only a broad flat band can be observed in the 0-0 correlations in the frequency range 
20 THz > z/ > 25 THz and in the Si-Si correlations an effective maximum appears around 
u = 20.5 THz. 

B. Density correlations 

Studying the density-density-correlation function in the g-i^-domain, i.e. the dynamic 
structure factor, 



Siq, ^) = j^ /_^ dt ( exp (z27rz/t) ^ exp (zq ■ (r,(t) - r,(0))) \ , (6) 

is of special interest because this quantity can be directly measured in neutron scattering 
experiments. S^q,^) is related to the longitudinal current correlation function Ji{q,h') by 
the simple equation 

which holds because of the continuity equation for particle number conservation [Q. Equa- 



tion (^ means that S{q, p) and Ji{q^ z/) contain the same information but features at lower 
frequencies are strongly enhanced in S{q, v) because of the factor 1/z/^. Moreover, S'(g, v) ex- 
hibits a quasielastic line around z/ = whereas J;(g, v) approaches zero for z/ ^ 0. Therefore, 
one has to investigate density fluctuations in order to understand the relationship between 
vibrational and relaxational dynamics, which is one of the goals of the present sectioni, since 
the latter is seen mainly at small v. 

In Fig. |lO| a S{q, v) is shown for several values of g at T = 2750 K. At the three lowest 
q values of our simulation, q = 0.13 A~^, 0.18 A~^ und 0.23 A~^, one sees essentially one 
peak which corresponds to the longitudinal acoustic excitations moving to higher frequencies 
with increasing q. At higher values of q the LA excitation is only visible as a shoulder until 
it reaches q = 1.7 A^^ at which it can be identified as a broad peak around z/ = 17 THz. 
The reason why this excitation is relatively hard to see is due to the fact that for q > 
0.23 A"^ a second peak is present in S{q, z/) which is located nearly independent of q around 
u = 1.7 THz. This peak is the so-called boson peak which is also seen experimentally for 
silica in Raman and neutron scattering 0,§]. From the figure we also see that to the left 
of the boson peak two additional peaks are present. The location of these two peaks is 
at u = 0.75 THz and u = 1.05 THz and we have checked that they are not due to bad 
statistics nor due to artifacts of the Fourier transformation. In order to discuss their origin 
we have plotted in Fig. p!0| b S{q, u) at q = 1.7 A'^ and Jt{q, z^) for the five lowest q values 
of our simulation q = 0.13 1"^ 0.18 A'^ 0.23 A'^ 0.26 A'^ and 0.29 A'^ Also included 
in Fig. |TD|b is the sum of these transverse current correlation functions Ji,sum (bold dashed 
line) which we have tried to scale onto S{q, v) by dividing it by 1.6 in order to compare the 
shape of this function with that of S{q, v). From the comparison of Ji,sum with the dynamic 
structure factor we can conclude that the main contribution to the low frequency part of 
the boson peak comes from the coupling to the TA modes at g = 0.13 A~^, 0.18 A~^, and 
0.23 A~^. The mechanism of how these modes couple to density fluctuations at higher g, 
e.g. at 1.7 A~^ in Fig. p!0|b, is due to elastic scattering since the energy of the scattered 
TA modes is conserved. We have seen before that the boson peak can only be observed for 
q > 0.23 A~^, which is exactly the region of q in which the LA and TA peaks in J{q,u) 
begin to overlap (see Fig. [l|a). That the transverse part is of special importance is plausible 
since the intensity of the TA peaks is a factor 6-8 higher than that of the LA peaks in the 
current correlation functions at fixed q (see Fig. p. Note that around z/ = 1.0 THz the band 
of acoustic modes is not dense in our simulation because of the finite size of the simulation 
box. For this reason one would expect that the intensity of the low frequency part of the 
boson peak is underestimated by our simulation. But as it is demonstrated by Fig. p!0| b this 
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property of the finite size system allows us to identify the influence of the low q TA modes 
on S{q, v) at much larger q and small v. 

One might argue that the experimental result, that the boson peak is observed in the 
Raman spectra of vitreous silica for q values around 10~^ A~^, is in contradiction to our 
simulation in which we find this peak in S{q, v) only for q > 0.23 A^^. But this is probably 
due to the fact that in Raman scattering a strong coupling of the light to the incoherent 
part of the density fluctuations is present. Indeed in our simulation the boson peak is also 
visible in the self part of the dynamic structure factor Ss{q, v) which is obtained from Eq. (^ 
for S{q, v) by taking into account only the terms with fc = / in the sum. As can be seen in 
Fig. p!0| c even at g = 0.13 A^^, the smallest accessible wave-vector, we observe the boson 
peak in S's((3', z/) and this quantity also exhibits the sharp peaks around v = 0.8 THz and 
z/ = 1.05 THz that stem from the TA modes. Furthermore, the shape of the different curves 
in this figure seems to be independent of q. If this is the case this means that Ss{q, i') can 
be factorized into a q dependent part and a purely frequency dependent part. Indeed this 



has been predicted recently in an analytic calculation by Gotze and Mayr [|T8[ for a hard 
sphere system within the mode coupling theory of the glass transition, and they found that 
the q dependent part is proportional to q^. In order to test whether this prediction holds 
we have plotted Ss{q,v)/q^ for silicon and oxygen in Figs. |TT]a and |TT]b, respectively. We 
recognize that the curves for 0.13 A~^ < 5' < 1.0 A~^ fall nicely onto one master curve in the 
whole frequency range 0.5 THz < z^ < 10 THz whereas at larger q small deviations from the 
master curve are visible around the location of the boson peak, z/ = 1.5 THz. To study the 
q dependence of Ss{q-, v) in more detail we show in Fig. |l^ a double logarithmic plot of this 
quantity at the frequencies 1.64 THz, 3.02 THz, 10.01 THz, and 30.02 THz. We see that fits 
with quadratic laws (bold lines in the figures) hold very well at least for q < 2.0 A~^. This 
means that the whole spectrum scales with g^ in this q range. Note that a similar behavior 
was also found in a simulation of ZnCl2 [P . 



From the harmonic approximation one would expect that S{q, u) scales with temperature. 
For this reason we have plotted in Fig. |13| S{q, z/)/T as a function of frequency at g = 1.7 A~^ 
for the temperatures T = 3760 K, 2750 K, 1670 K, 1050 K, and 300 K. We recognize that 
the curves for the different temperatures fall roughly on top of each other. This means that 
also in the region of the boson peak our silica model is quite harmonic even at the relatively 
high temperature T = 3760 K. Note that even at this temperature the boson peak feature 
is present in that a shoulder is formed for frequencies above 0.5 THz. 

If it is indeed true that TA modes with q < 0.2 A~^ couple to density fluctuations at 
higher q giving rise to certain features in the boson peak, such as the additional peak at 
z/ = 0.8 THz, then these features should be absent if the system size is so small that it does 
not allow the propagation of TA modes with q < 0.2 A~^. In order to check this prediction 
we have calculated Ss{q, z/) at T = 3760 K for the system sizes N = 336 and A^ = 1002 in 
addition to A^ = 8016. As the same density as for A^ = 8016 is used, pm = 2.37 g/cm^, 
the sizes of the simulation boxes are L = 16.80 A and L = 24.18 A for A^ = 336 and 
A^ = 1002, respectively. Thus the smallest q values are q = 0.37 A^^ and q = 0.26 A^^. In 
Fig. 0a we show the obtained Ss{q, z/) at g = 0.37 A~^, 1.7 A~^, and 4.75 A~^ for the three 
system sizes. Whereas the curves for the different system sizes coincide for frequencies that 
are larger than a weakly A^ dependent frequency z/cut(A^), for z/ < z/cut(A^) the magnitude 
of 5s(g, z/) decreases with decreasing A^. Independent of g we read off z/cut ~ 1-7 THz for 
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A^ = 336 and z/cut ~ 1.2 THz for A^ = 1002. Both frequencies are marked by vertical lines 
in Fig. p!4|a. fcut(-^) is indeed just slightly below the frequency of the transverse excitation 
corresponding to the lowest q value determined by the size of the simulation box. These 
frequencies are aX v = 1.85 THz for A^ = 336 and at u = 1.35 THz for A^ = 1002. Thus this 
is evidence that the missing of the TA modes with q < 0.2 A~^ causes the finite size effects 
in the small systems. 

Due to the sum rule J du Ss{q,i^) = 1, the loss of intensity in the boson peak below 
^'cut(^) has to be "reshuffled" to smaller frequencies leading to a broadening and an increase 
of the intensity of the quasielastic line around z/ = 0. Since the quasielastic line is outside 
the frequency resolution of our Fourier transformation the consequences in the change of 
the quasielastic line can be observed better in the Fourier transform of S's(g, z/), i.e. the 
incoherent intermediate scattering function Fs{q,t) which is defined as 

AT foo 

Fs{q, t) = ^ du exp (27rz/t) S,{q, u) a G [Si, O] . (8) 

iV J— oo 

For oxygen this quantity is shown in Fig. |l^ b at g = 1.7 A^^ for different system sizes. (We 
mention that we have also included a curve for A^ = 3006 in the figure. The box length in 
this case is L = 34.89 A corresponding to the lowest q value q = 0.18 A~^.) We recognize 
from Fig. |T^ that with decreasing system size the height of the plateau, which is the Lamb- 
Mossbauer factor, increases and the a-relaxation time shifts to longer times. Furthermore, 
the scattering functions for the small systems show a pronounced oscillation for t > 0.2 ps. 
This can be simply understood from the z^— dependence of S's(g, z/): For A^ = 8016 this 
quantity shows a shoulder around 1.0 THz which corresponds to the monotonous decay of 
Fs{q,t). In the small systems there is a peak in Ss{q, z/) around z/cut(^) which corresponds 
to the oscillations in Fs{q,t) with a period l/z/cut(A^). From the fact that the band of the 
transverse acoustic modes is not dense for the region of small q (see Fig. p!0|c) we expect also 



for A^ = 8016 that the finite size effects are present. But, since in Fig. |Tjb the differences of 
the curves for A^ = 3006 and A^ = 8016 are small, we can conclude that the finite size effects 
play not an important role for A^ = 8016, at least for T = 3760 K. In order to investigate 
the influence of these flnite size effects on the a— relaxation we deflne the a-relaxation time 
Ta as the time at which Fs{q, t) has decayed to 1/e. Ta increases from A^ = 8016 to A^ = 336 
by about 40 %. However, in the a-relaxation regime the shape of Fs{q,t) does not seem to 
depend on the system size, as can be seen in the inset of Fig. |1^ where we have plotted 
Fs{q,t) versus the scaled time t/r^. We see that in the a-relaxation regime the curves for 
the different system sizes fall onto one master curve. This holds also for Fs{q,t) for the 
silicon atoms. 

Of course, the flnite size effects are also important in the total dynamic structure factor. 



Figure p!5| a shows S{q, u) for the two system sizes A^ = 8016 and A^ = 336 at T = 2750 K 
and the three q values q = 0.9 A~^, 1.7 A~^ and 4.75 A^^. We can again identify a cut-off 
frequency around 1.7 THz below which there is a loss of intensity in S{q, v) for A^ = 336. 
Note that the sharp peaks at z/ = 0.75 THz and z/ = 1.05 THz are again not present in the 
small system. Moreover, we recognize that the relative intensity loss in the small system 
compared to the large system depends on q. In order to quantify this q dependence we deflne 
the ratio 

R{(l,^)-=^ 7 r-1 (9) 
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which is zero if the dynamic structure factor coincides for the two system sizes. Figure [15[ b 
shows R{q, v) for several values of q. Its behavior underlines what we have said before that 
the low frequency part of the boson peak is mainly due to the elastic scattering of the two 
TA modes with the frequencies v = 0.75 THz and u = 1.05 THz corresponding to the lowest 
q values of our simulation for the system size A^ = 8016. Obviously, the amplitudes of the 
peaks in R{q, v) do not change monotonously as a function of q. In order to investigate 
this q dependence i?(g, v) is plotted in Fig. |T^ for the frequencies v = 0.75 THz and 
z/ = 1.05 THz as a function of q. R{q, v) shows pronounced maxima at g = 1.6 A~^, i.e. in 
the vicinity of the location of the first sharp diffraction peak in the static structure factor, 
and at g = 2.8 A~^, which is the location of the second peak in S{q, u). Thus, the structural 
disorder on intermediate length scales, i.e. the length scale introduced by two connected 
Si04-tetrahedra, is most relevant for the scattering of the TA modes with q < 0.2 A~^. Also 
included is Rs{q, v) for the two frequencies which is obtained by putting in 5*8 (g, z/) instead 
of S{q,u) into the definition (^). The incoherent part Rs^q^u) decreases monotonously 
with q which is plausible since the finite size effects should vanish at very large values of q 
corresponding to small length scales. 

IV. SUMMARY AND CONCLUSIONS 

We have done molecular dynamics simulations in order to investigate the high frequency 
dynamics of amorphous silica. The results which we have presented in this paper have been 
for the fully equihbrated liquid and the glass state. The frequency range we have studied is 
0.5 THz < V < AQ THz for wave-vectors with magnitude g > 0.13 A~^ (limited by the size 
of the simulation box). 

In a first step we have discussed the properties of the longitudinal and transverse current 
correlation functions. At low frequencies we have identified propagating longitudinal acoustic 
(LA) and transverse acoustic (TA) modes the maxima of which move to higher frequencies 
with increasing g. The amplitude of the TA peaks is a factor 6-8 larger than that of the LA 
peaks at a fixed value of g which is a first indication for the importance of the transverse 
dynamics in silica. Whereas the LA peak is separated quite well from the TA peak on the 
frequency axis at g = 0.13 A~^, both peaks begin to overlap at higher g. The g region at 
which the LA and TA peaks begin to overlap significantly can be seen as a crossover region 
from a regime where the longitudinal and transverse modes exhibit only a weak interaction 
to a regime where a strong interaction between different modes is present. One important 
sign of this is that the qualitative shape of the spectra starts to change gradually around 
g = 0.6 A~^: The LA peaks are still well-defined, but they are now sit on top of a flat 
background. The acoustic band in Jf (g, z/) shows a similar behavior in that it evolves into a 
broad plateau from about 3 to 11 THz. The observation that the acoustic modes are located 
on top of a fiat background for intermediate values of g has recently been found by Gotze 
and Mayr |]18| as an essential result in an analytic calculation of the spectra of hard sphere 



systems in their glass state within the mode-coupling theory of the glass transition. Within 
their theory Gotze and Mayr have explained the existence of the fiat background spectrum 
by the presence of inelastic phonon scattering where a mode decays into two modes due to 
anharmonicity. In the same sense we can understand the fiat background spectrum in silica 
as the manifestation of multiphonon excitations. 
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By reading off the peak maxima in Ji{q, v) and Jt{<ii v) as a function of q one obtains 
dispersion like functions Viiq) and Vtiq) for the longitudinal and transverse part, respectively. 
z//(g) shows an approximately linear behavior for wave-vectors up to 0.3 A~^. On the other 
hand the full width at half maximum Ti{q) of the LA peaks is well decribed by a quadratic 
law for q < 0.5 A^^. This means that to a good approximation the system behaves like an 
isotropic elastic medium up to g ~ 0.3 A^^ with respect to the longitudinal sound modes. 
Furthermore, vi{q) exhibits a quasi Brillouin zone at gm/2 where q^ is the location of the 
second sharp diffraction peak in the static structure factor correponding to length scales 
of intratetrahedral distances. Also Vtiq) shows approximately a linear behavior at small q. 
But the data for Tt{q) cannot be described by a quadratic law. Instead, Tt{q) is fitted well 
with a q^-^ law from which we conclude that the TA excitations are stronger damped than 
expected from an isotropic elastic medium. For q > 0.8 A~^ ^tiq) becomes flat, and in the 
same range the TA excitations become strongly overdamped in that they reach a loffe-Regel 
limit, i.e. Tt{q) is of the order of z/j(g). We emphasize that there are no pecularities in the 
qualitative behavior of t^iiq) and I'tiq) for our silica model since these functions look very 
similar for simple liquids [Q. 

From the two lowest q values of our simulation, q = 0.13 A~^ and 0.18 A~^, we have 
determined the apparent high frequency sound velocities for the different temperatures and 



find that they reproduce the light scattering data by Vo-Tanh et al. ||38| very well. Thus 
this is another example that the BKS model is able to reproduce reliably the experimental 
data of amorphous silica. 

Most of the results we have summarized up to now for the acoustic band are for the 
silica melt at T = 2750 K. But we have seen that in the temperature range 300 K < 
T < 3760 K the spectra scale to a good approximation linearly with temperature for 
frequencies v < 2Q THz, a behavior which is expected if the harmonic approximation is 
valid. The linear scaling is not valid for the complex optical bands above 20 THz. For 
these excitations the spectrum becomes more pronounced with decreasing temperature. If 
one plots the current correlation functions divided by temperature one recognizes that the 
peak maxima corresponding to the intratetrahedral stretching modes L02 and T02 move 
to higher frequencies and also their amplitude increases with decreasing temperature. This 
is just the result of the fact that the particles are more localized the lower the temperature 
is, which causes anharmonicities to dissapear. Nevertheless, due to the disorder, especially 
for the L02 and T02 modes, a strong mixing occurs between the longitudinal and the 
transverse part for q > 0.3 A~^ which can be clearly identified at T = 300 K. So, these 
wave-vectors are no longer good quantum numbers, and if one treats disordered structures 
within the harmonic approximation one has to keep in mind that this approximation only 
describes the "mixing contributions" in an effective way. 

In a second step we have discussed density fluctuations by means of the dynamic structure 
factor S{q^ u). For q > 0.23 A~^ this quantity exhibits a boson peak which is located nearly 
q independent around z^bp = 1-7 THz at T = 2750 K. The boson peak excitations coexist 
with the LA modes since the latter is visible also at frequencies above z/Bp. Since the boson 
peak has a much larger intensity than the LA peak, e.g. a factor of 7-8 for g = 1.0 A~^, the 
LA excitations are visible only as a shoulder in S{q, v). Only if the LA peak has moved to 
high frequencies, e.g. to 17 THz at g = 1.7 A~^, one observes this peak as a second peak in 
addition to the boson peak in the dynamic structure factor. In the low frequency part of the 
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boson peak two sharp peaks are present aX v = 0.75 THz and v = 1.05 THz which are due 
to the elastic scattering of the TA modes with q = 0.13 A~^ and q = 0.18 A~^, respectively. 
We will discuss them in more detail below. Also the dynamic structure factor S{q, z/) scales 
for frequencies around u^p roughly with temperature in the range 3760 K > T > 300 K. 
Note that at T = 3760 K the boson peak feature can be only seen as a shoulder which 
grows out of the quasielastic line. The fact that the boson peak can be seen even at such 
high temperatures supports the view of Ref. [^ that this feature becomes visible as soon 
as temperature is around T^, which for our system is around 3330 K [p^] . 

For wave- vectors up to about 1-2 A~^ the self part of the dynamic structure factor 
Ss{q, v) exhibits a factorization into a frequency dependent and wave- vector dependent part 
whereby the latter is proportional to q^. Apart from the fact that this property of Ss{q-, v) 
has also been found in a MD simulation of ZnCl2 ^1|] it is remarkable that this result comes 



out of the mode coupling calculation of Goetze and Mayr [|18|. So this is another important 
feature which is reproduced by this theoretical approach. 

In order to get more insight into the boson peak feature we have done simulations also 
for smaller system sizes than the normally used system with A^ = 8016 particles. We have 
found strong finite size effects in the low frequency part of the boson peak which can be 
characterized by a frequency z/cut(A^) below which there is a lack of intensity in the dynamic 
structure factor. The frequency Vc\it{N) decreases with increasing system size N and is 
essentially independent of q. The reason for these finite size effects is due to the absence of 
the TA excitations with q < 0.2 A^^ in the small systems since the smallest q value of our 
simulations with A^ = 1002 and A^ = 336 particles are 0.26 A~^ and 0.37 A~^, respectively. 
In the time domain, i.e. in the incoherent intermediate scattering function Fs{q,t), the finite 
size effects cause an increase of the Lamb-Mossbauer factor and of the a— relaxation time. 
This is a consequence of the sum rule J du 5*8 (g,//) = 1 since the missing of intensity for 
^ < t'cut(A^) has to be "reshuffled" to smaller frequencies. Because of the abrupt decrease of 
Ss{q, v) below t'cut(A^) for small A^ we observe quite pronounced oscillations for t > 0.2 ps in 
Fs{q, t) whereby the period of these oscillations is approximately equal to z/cut(A^)- Note that 



a similar behavior was also found in a MD simulation by Lewis and Wahnstrom [42] for a 



model of orthoterphenyl in which the interactions between the molecules are described by a 
Lennard- Jones potential. These authors have suggested that a disturbance that propagates 
through the system will leave and reenter the box due to the periodic boundary conditions 
after a time L/c, where L is the size of the box and c is the typical velocity of the sound 
wave. This mechanism then produces an echo, i.e. an additional signal which produces the 
slowed down decay of the correlation functions like Fs{q,t). We have also suggested this 
explanation recently for silica [^], but we think now that this explanation for the finite 
size effects is not the correct one. Instead, the general argument is that in small enough 
systems (with a smallest wave-vector with magnitude gs) parts of the vibrational spectrum 
are missing below a frequency z/^ut ( A^) because of the coupling to wave-vectors with q < qg 
occuring in an infinite system. In a Lennard- Jones system such a coupling is reflected in 
the flat background spectrum which was shown to be present by Mazzacurati et al. |Q. In 
the case of silica there is in addition the coupling which arises from the elastic scattering of 
transverse TA modes with small q by the structural disorder. 

One may speculate that the strength of the boson peak in silica is due the strong coupling 
of the TA exciations to the longitudinal part. The stiffness of the tetrahedral Si02 network 
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introduces strong restoring forces which allow the propagation of shear waves with a large 
amplitude even at relatively high temperatures. We have shown that the strongest scattering 
of the TA modes is at g ~ 1.6 A~^. This is in agreement with suggestions in the literature 
that the boson peak is caused by the interactions of sound modes with coupled rotations of 
several tetrahedra P,|35|. A possible mechanism of this interaction would be that the coupled 
rotations of the tetrahedra enable the change of the polarization of transverse acoustic modes, 
so that they contribute to the density fluctuation spectrum, i.e. constituting at least part of 
the boson peak. 
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FIG. 1. Longitudinal and transverse current correlation functions (filled and open symbols, 
respectively) for the oxygen-oxygen correlations at the temperature T = 2750 K. a) The peaks 
moving to higher frequencies correspond to q = 0.13 A~^, q = 0.18 A^^, q = 0.23 A~^, and 
q = 0.26 A~^ for the longitudinal and transverse functions, respectively. Note that the curves for 
Ji{q, u) are multiplied by a factor of 4. h) q = 1.0 A^^. The solid line corresponds to the sum of 
Ji and Jt- 
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FIG. 2. Ji{q,v) and Jt{q.,y) for various wave-vectors at T = 2750 K, a) Ji for 0-0, b) Ji for 

Si-0, c) Ji for Si-Si, d) Jt for 0-0, e) Jt for Si-0, f) Jj for Si-Si. 
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FIG. 3. Peak maximuiii position fi^tio) for the LA and TA modes for the Si-Si correlations 
(filled symbols) and the 0-0 correlations (open symbols) at a) T = 2750 K and b) T = 300 K. 
The bold lines are fits with linear laws Va = Co-g/ (27r) for which the corresponding values for the 
sound velocities Ca are given in the figure. 
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FIG. 4. The temperature dependence of the sound velocities ci^t which are determined from 
lyi^til) at g = 0.13 A~^ (open circles) and q = 0.18 A^^ (filled squares). Also included are the 
experimental data of Vo-Tanh et al. |g^ which we have multiplied with the factor (2.2/2.37)'^'^ in 
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FIG. 5. a) Ji{q,i^)/Ji^ina.x versus ly — t'/(g) at T = 2750 K for the 0-0 correlations in the q 
range 0.13 A"! < q < 2.0 k''^. b) Jt{q,iy)/Jt,mi,K versus u - ut{q) at T = 2750 K for the 0-0 
correlations in the q range 0.13 A~^ < ^7 < 0.5 A~^. 
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FIG. 6. Linewidth Ti^tiQ) obtained at T = 2750 K. The bold lines are power law fits with 
exponents 2 and 2.5 for Ti{q) and Tf{q), respectively. Also included are lyiiq) (bold dashed line) 
and t't(Q') (bold dotted line). 
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FIG. 7. The current correlation functions scaled with temperature for q = 0.26 A~^, a) 
Ji{q,u)/T for the 0-0 correlations, b) Jt{q,i')/T for the 0-0 correlations, c) Ji{q,v)/T for the 
Si-O correlations, d) Jt{q,i')/T for the Si-O correlations, e) Ji{q,v)/T for the Si-Si correlations, 
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FIG. 8. The same as in Fig. but now for q = 1.7 K 
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FIG. 9. Peak maximum position vi^t{(l) of the excitations in the optical band for the Si-Si 
correlations (filled symbols) and the 0-0 correlations (open symbols) at a) T = 300 K and b) 
T = 2750 K. See text for the explanation of the dashed and the solid line in a). 
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FIG. 10. a) Dynamic structure factor S{q, v) for the 0-0 correlations at T = 2750 K for several 
values of q. b) S{q,v) at q = 1.7 A~^ (bold solid line) and Jt{q,i') for the five lowest q values 
of our simulation. The bold dashed line shows the sum of the latter current correlation functions 
divided by 1.6. c) Self part of the dynamic structure factor Ss{q, i^) for q = 0.13 A~^, q = 0.23 A~^, 
q = 0.32 A-\ and q = 0.41 A"!. 
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FIG. 11. Self part of the dynamic structure factor divided by q^ for q = 0.13 A ^, 0.6 A ^, 
1.0 A^^, 1.7 A~^, and 2.75 A~^ for a) Si and b) O. The vertical lines are at the frequencies 



u = 1.64 THz and 3.02 THz at which the q dependence of S{q, u) is shown in Fig. 12. 
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FIG. 12. q dependence of Ss{q,y) for v = 1.64 THz, 3.02 THz, 10.01 THz, and 30.02 THz, for 
a) Si, and b) O. The bold straight hnes are fits with q^ laws. 
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FIG. 13. S{q,v)/T for the 0-0 correlations at g = 1.7 A ^ for the different temperatures. 
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FIG. 14. a) The self part of the dynamic structure factor for oxygen for the different system 
sizes N = 8016, 1002,336 at the wave-vectors q = 0.37 A'^ 1.7 A'^, and 4.75 A'^, T = 3760 K. 
See text for the explanation of the vertical lines, b) Incoherent intermediate scattering function for 
oxygen for the system sizes N = 8016, 3006, 1002, and 336 at the wave-vector q = 1.7 A^^ and the 
temperature T = 3760 K. The inset shows the same data (without the N = 3006 curve), plotted 
versus the scaled time t/ra- 
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FIG. 15. a) The dynamic structure factor for the 0-0 correlations for the system sizes 
A^ = 8016 and N = 336 at the wave-vectors q = 0.9 1"^ 1.7 A'^, and 4.75 A'^ and the 
temperature T = 2750 K. b) R{q, f) (definition see text) for several values of g as a function of fre- 
quency at T = 2750 K. c) R{q, v) for v = 0.75 THz (filled circles) and v = 1.05 THz (open squares) 
as a function of g at T = 2750 K. Also included is Rs{q, v) at the same frequencies obtained from 
Ssiq, v) (bold lines). 
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